home *** CD-ROM | disk | FTP | other *** search
/ Amiga Format CD 46 / Amiga Format CD46 (1999-10-20)(Future Publishing)(GB)[!][issue 1999-12].iso / -in_the_mag- / reader_requests / scilab / tests / dassldasrt.dia.ref < prev    next >
Text File  |  1999-09-16  |  14KB  |  476 lines
  1.  
  2. //DASSL
  3.  
  4. // PROBLEM 1..   LINEAR DIFFERENTIAL/ALGEBRAIC SYSTEM
  5.  
  6. //
  7.  
  8. //X1DOT + 10.0*X1 = 0
  9.  
  10. //X1 + X2 = 1
  11.  
  12. //X1(0) = 1.0, X2(0) = 0.0
  13.  
  14. //
  15.  
  16. t=1:10;t0=0;y0=[1;0];y0d=[-10;0];
  17.  
  18. info=list([],0,[],[],[],0,0);
  19.  
  20. //    Calling Scilab functions
  21.  
  22. deff('[r,ires]=dres1(t,y,ydot)','r=[ydot(1)+10*y(1);y(2)+y(1)-1];ires=0')
  23.  
  24. comp(dres1);
  25.  
  26. deff('pd=djac1(t,y,ydot,cj)','pd=[cj+10,0;1,1]')
  27.  
  28. comp(djac1);
  29.  
  30. //   scilab function, without jacobian
  31.  
  32. yy0=dassl([y0,y0d],t0,t,dres1,info);
  33.  
  34. //   scilab functions, with jacobian
  35.  
  36. yy1=dassl([y0,y0d],t0,t,dres1,djac1,info);
  37.  
  38. // fortran routine, without jocabian
  39.  
  40. yy2=dassl([y0,y0d],t0,t,'dres1',info);   //=yy0
  41.  
  42. if norm(yy2-yy0,1)>1D-5 then bugmes();quit;end
  43.  
  44. // fortran routines, with jacobian
  45.  
  46. yy3=dassl([y0,y0d],t0,t,'dres1','djac1',info);  //=yy1
  47.  
  48. if norm(yy3-yy1,1)>1D-5 then bugmes();quit;end
  49.  
  50. yy3bis=dassl([y0,y0d],t0,t,'dres1',djac1,info);
  51.  
  52. // call fortran dres1 and scilab's djac1
  53.  
  54. yy3ter=dassl([y0,y0d],t0,t,dres1,'djac1',info);
  55.  
  56. //
  57.  
  58. // with specific atol and rtol parameters
  59.  
  60. atol=1.d-6;rtol=0;
  61.  
  62. yy4=dassl([y0,y0d],t0,t,atol,rtol,dres1,info);
  63.  
  64. yy5=dassl([y0,y0d],t0,t,atol,rtol,'dres1',info); //=yy4
  65.  
  66. if norm(yy5-yy4,1)>1D-9 then bugmes();quit;end
  67.  
  68. yy6=dassl([y0,y0d],t0,t,atol,rtol,dres1,djac1,info);
  69.  
  70. yy7=dassl([y0,y0d],t0,t,atol,rtol,'dres1','djac1',info); //==yy6
  71.  
  72. if norm(yy7-yy6,1)>1D-12 then bugmes();quit;end
  73.  
  74. //
  75.  
  76. //   Testing E xdot - A x=0
  77.  
  78. //   x(0)=x0;   xdot(0)=xd0
  79.  
  80. rand('seed',0);
  81.  
  82. nx=5;
  83.  
  84. E=rand(nx,1)*rand(1,nx);A=rand(nx,nx);
  85.  
  86. //         Index 1
  87.  
  88. [Si,Pi,Di,o]=penlaur(E,A);pp=Si*E;[q,m]=fullrf(pp);x0=q(:,1);x0d=pinv(E)*A*x0;
  89.     rank A^k    rcond
  90.         1.      0.100D+01
  91.  
  92. deff('[r,ires]=g(t,x,xdot)','r=E*xdot-A*x;ires=0');comp(g);
  93.  
  94. t=[1,2,3];t0=0;info=list([],0,[],[],[],0,0);
  95.  
  96. x=dassl([x0,x0d],t0,t,g,info);x1=x(2:nx+1,:);
  97.  
  98. if norm(pp*x1-x1,1)>1.d-5 then bugmes();quit;end
  99.  
  100. //x(4) goes through 1 at  t=1.5409711;
  101.  
  102. t=1.5409711;ww=dassl([x0,x0d],t0,t,g,info);
  103.  
  104. if abs(ww(5)-1)>0.001 then bugmes();quit;end
  105.  
  106. deff('[rt]=surface(t,y,yd)','rt=y(4)-1');nbsurf=1;
  107.  
  108. [yyy,nnn]=dasrt([x0,x0d],t0,t,g,nbsurf,surface,info);
  109.  
  110. deff('pd=j(t,y,ydot,cj)','pd=cj*D-A');comp(j);
  111.  
  112. x=dassl([x0,x0d],t0,t,g,j,info);x2=x(2:nx+1,1);
  113.  
  114. if norm(x2-ww(2:nx+1,1),1)>0.0001 then bugmes();quit;end
  115.  
  116. [yyy1,nnn]=dasrt([x0,x0d],t0,t,g,j,nbsurf,surface,info);
  117.  
  118. //x0d is not known:
  119.  
  120. x0d=ones(x0);info(7)=1;
  121.  
  122. x=dassl([x0,x0d],t0,t,g,info);
  123.  
  124. xn=dassl([x0,x0d],t0,t,g,j,info);
  125.  
  126. if norm(x-xn,1)>0.00001 then bugmes();quit;end
  127.  
  128.  
  129.  
  130. //PROBLEM 2..
  131.  
  132.  
  133. info=list([],0,[],[],[],0,0);
  134.  
  135. y0=zeros(25,1);y0(1)=1;
  136.  
  137. delta=0*y0;
  138.  
  139. //link('dres2.o','dres2');
  140.  
  141. //y0d=fort('dres2',0,1,'d',y0,2,'d',delta,3,'d',0,5,'i',0,6,'d',0,7,'d','out',[25,1],4,'d');
  142.  
  143. y0d=zeros(y0);y0d(1)=-2;y0d(2)=1;y0d(6)=1;
  144.  
  145. t0=0;t=[0.01,0.1,1,10,100];
  146.  
  147. rtol=0;atol=1.d-6;
  148.  
  149. y=dassl([y0,y0d],t0,t,atol,rtol,'dres2',info);
  150.  
  151.  
  152. //                 DASRT
  153.  
  154. //
  155.  
  156. //C-----------------------------------------------------------------------
  157.  
  158. //C First problem.
  159.  
  160. //C The initial value problem is..
  161.  
  162. //C   DY/DT = ((2*LOG(Y) + 8)/T - 5)*Y,  Y(1) = 1,  10.LE. T0.LE. 6
  163.  
  164. //C The solution is  Y(T) = EXP(-T**2 + 5*T - 4), YPRIME(1) = 3
  165.  
  166. //C The two root functions are..
  167.  
  168. //C   G1 = ((2*LOG(Y)+8)/T - 5)*Y (= DY/DT)  (with root at T = 2.5),
  169.  
  170. //C   G2 = LOG(Y) - 2.2491  (with roots at T = 2.47 and 2.53)
  171.  
  172. //C-----------------------------------------------------------------------
  173.  
  174. y0=1;t=2:6;t0=1;y0d=3;
  175.  
  176. info=list([],0,[],[],[],0,0);
  177.  
  178. atol=1.d-6;rtol=0;ng=2;
  179.  
  180. [yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,'res1',ng,'gr1',info);
  181.  
  182. if abs(nn(1)-2.47)>0.001 then bugmes();quit;end
  183.  
  184. y0=yy(2,2);y0d=yy(3,2);t0=nn(1);t=[3,4,5,6];
  185.  
  186. [yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,'res1',ng,'gr1',info);
  187.  
  188. if abs(nn(1)-2.5)>0.001 then bugmes();quit;end
  189.  
  190. y0=yy(2,1);y0d=yy(3,1);t0=nn(1);t=[3,4,5,6];
  191.  
  192. [yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,'res1',ng,'gr1',info);
  193.  
  194. if abs(nn(1)-2.53)>0.001 then bugmes();quit;end
  195.  
  196.  
  197. deff('[delta,ires]=res1(t,y,ydot)','ires=0;delta=ydot-((2*log(y)+8)/t-5)*y')
  198.  
  199. deff('[rts]=gr1(t,y,yd)','rts=[((2*log(y)+8)/t-5)*y;log(y)-2.2491]')
  200.  
  201. comp(res1);comp(gr1);
  202.  
  203.  
  204. y0=1;t=2:6;t0=1;y0d=3;
  205.  
  206. info=list([],0,[],[],[],0,0);
  207.  
  208. atol=1.d-6;rtol=0;ng=2;
  209.  
  210. [yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,res1,ng,gr1,info);
  211.  
  212. if abs(nn(1)-2.47)>0.001 then bugmes();quit;end
  213.  
  214. y0=yy(2,2);y0d=yy(3,2);t0=nn(1);t=[3,4,5,6];
  215.  
  216. [yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,res1,ng,gr1,info);
  217.  
  218. if abs(nn(1)-2.5)>0.001 then bugmes();quit;end
  219.  
  220. y0=yy(2,1);y0d=yy(3,1);t0=nn(1);t=[3,4,5,6];
  221.  
  222. [yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,res1,ng,gr1,info);
  223.  
  224. if abs(nn(1)-2.53)>0.001 then bugmes();quit;end
  225.  
  226.  
  227. //C
  228.  
  229. //C-----------------------------------------------------------------------
  230.  
  231. //C Second problem (Van Der Pol oscillator).
  232.  
  233. //C The initial value problem is..
  234.  
  235. //C   DY1/DT = Y2,  DY2/DT = 100*(1 - Y1**2)*Y2 - Y1,
  236.  
  237. //C   Y1(0) = 2,  Y2(0) = 0,  00.LE. T0.LE. 200
  238.  
  239. //C   Y1PRIME(0) = 0, Y2PRIME(0) = -2
  240.  
  241. //C The root function is  G = Y1.
  242.  
  243. //C An analytic solution is not known, but the zeros of Y1 are known
  244.  
  245. //C to 15 figures for purposes of checking the accuracy.
  246.  
  247. //C-----------------------------------------------------------------------
  248.  
  249. rtol=[1.d-6;1.d-6];atol=[1.d-6;1.d-4];
  250.  
  251. t0=0;y0=[2;0];y0d=[0;-2];t=[20:20:200];ng=1;
  252.  
  253. info=list([],0,[],[],[],0,0);
  254.  
  255. [yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,'res2','jac2',ng,'gr2',info);
  256.  
  257. if abs(nn(1)-81.163512)>0.001 then bugmes();quit;end
  258.  
  259.  
  260. deff('[delta,ires]=res2(t,y,ydot)',...
  261. 'ires=0;y1=y(1),y2=y(2),delta=[ydot-[y2;100*(1-y1*y1)*y2-y1]]')
  262.  
  263. [yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,res2,'jac2',ng,'gr2',info);
  264.  
  265. deff('J=jac2(t,y,ydot,c)','y1=y(1);y2=y(2);J=[c,-1;200*y1*y2+1,c-100*(1-y1*y1)]')
  266.  
  267. [yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,res2,jac2,ng,'gr2',info);
  268.  
  269. deff('s=gr2(t,y,yd)','s=y(1)')
  270.  
  271. [yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,res2,jac2,ng,gr2,info);
  272.  
  273.  
  274. //           Hot Restart
  275.  
  276.  
  277. [yy,nn,hot]=dasrt([y0,y0d],t0,t,atol,rtol,'res2','jac2',ng,'gr2',info);
  278.  
  279. t01=nn(1);t=100:20:200;[pp,qq]=size(yy);y01=yy(2:3,qq);y0d1=yy(3:4,qq);
  280.  
  281. [yy,nn,hot]=dasrt([y01,y0d1],t01,t,atol,rtol,'res2','jac2',ng,'gr2',info,hot);
  282.  
  283. if abs(nn(1)-162.57763)>0.001 then bugmes();quit;end
  284.  
  285.  
  286. //Test of Dynamic link (Require f77!)
  287.  
  288. //         1 making the routines
  289.  
  290. res22=[...
  291. '      SUBROUTINE RES22(T,Y,YDOT,DELTA,IRES,RPAR,IPAR)';
  292. '      IMPLICIT DOUBLE PRECISION (A-H,O-Z)';
  293. '      INTEGER NEQ';
  294. '      DIMENSION Y(*), YDOT(*), DELTA(*)';
  295. '      NEQ=2';
  296. '      CALL F2(NEQ,T,Y,DELTA)';
  297. '      DO 10 I = 1,NEQ';
  298. '         DELTA(I) = YDOT(I) - DELTA(I)';
  299. ' 10   CONTINUE';
  300. '      RETURN';
  301. '      END';
  302. '      SUBROUTINE F2 (NEQ, T, Y, YDOT)';
  303. '      IMPLICIT DOUBLE PRECISION (A-H,O-Z)';
  304. '      INTEGER NEQ';
  305. '      DOUBLE PRECISION T, Y, YDOT';
  306. '      DIMENSION Y(*), YDOT(*)';
  307. '      YDOT(1) = Y(2)';
  308. '      YDOT(2) = 100.0D0*(1.0D0 - Y(1)*Y(1))*Y(2) - Y(1)';
  309. '      RETURN';
  310. '      END';]
  311.  res22  =
  312.  
  313. !      SUBROUTINE RES22(T,Y,YDOT,DELTA,IRES,RPAR,IPAR)    !
  314. !                                                         !
  315. !      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                !
  316. !                                                         !
  317. !      INTEGER NEQ                                        !
  318. !                                                         !
  319. !      DIMENSION Y(*), YDOT(*), DELTA(*)                  !
  320. !                                                         !
  321. !      NEQ=2                                              !
  322. !                                                         !
  323. !      CALL F2(NEQ,T,Y,DELTA)                             !
  324. !                                                         !
  325. !      DO 10 I = 1,NEQ                                    !
  326. !                                                         !
  327. !         DELTA(I) = YDOT(I) - DELTA(I)                   !
  328. !                                                         !
  329. ! 10   CONTINUE                                           !
  330. !                                                         !
  331. !      RETURN                                             !
  332. !                                                         !
  333. !      END                                                !
  334. !                                                         !
  335. !      SUBROUTINE F2 (NEQ, T, Y, YDOT)                    !
  336. !                                                         !
  337. !      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                !
  338. !                                                         !
  339. !      INTEGER NEQ                                        !
  340. !                                                         !
  341. !      DOUBLE PRECISION T, Y, YDOT                        !
  342. !                                                         !
  343. !      DIMENSION Y(*), YDOT(*)                            !
  344. !                                                         !
  345. !      YDOT(1) = Y(2)                                     !
  346. !                                                         !
  347. !      YDOT(2) = 100.0D0*(1.0D0 - Y(1)*Y(1))*Y(2) - Y(1)  !
  348. !                                                         !
  349. !      RETURN                                             !
  350. !                                                         !
  351. !      END                                                !
  352.  
  353.  
  354. jac22=[...
  355. '      SUBROUTINE JAC22 (T, Y, ydot, PD, CJ, RPAR, IPAR)';
  356. ' ';
  357. '      IMPLICIT DOUBLE PRECISION (A-H,O-Z)';
  358. '      INTEGER  NROWPD';
  359. '      DOUBLE PRECISION T, Y, PD';
  360. '      PARAMETER (NROWPD=2)';
  361. '      DIMENSION Y(2), PD(NROWPD,2)';
  362. '      PD(1,1) = 0.0D0';
  363. '      PD(1,2) = 1.0D0';
  364. '      PD(2,1) = -200.0D0*Y(1)*Y(2) - 1.0D0';
  365. '      PD(2,2) = 100.0D0*(1.0D0 - Y(1)*Y(1))';
  366. '      PD(1,1) = CJ - PD(1,1)';
  367. '      PD(1,2) =    - PD(1,2)';
  368. '      PD(2,1) =    - PD(2,1)';
  369. '      PD(2,2) = CJ - PD(2,2)';
  370. '      RETURN';
  371. '      END';]
  372.  jac22  =
  373.  
  374. !      SUBROUTINE JAC22 (T, Y, ydot, PD, CJ, RPAR, IPAR)  !
  375. !                                                         !
  376. !                                                         !
  377. !                                                         !
  378. !      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                !
  379. !                                                         !
  380. !      INTEGER  NROWPD                                    !
  381. !                                                         !
  382. !      DOUBLE PRECISION T, Y, PD                          !
  383. !                                                         !
  384. !      PARAMETER (NROWPD=2)                               !
  385. !                                                         !
  386. !      DIMENSION Y(2), PD(NROWPD,2)                       !
  387. !                                                         !
  388. !      PD(1,1) = 0.0D0                                    !
  389. !                                                         !
  390. !      PD(1,2) = 1.0D0                                    !
  391. !                                                         !
  392. !      PD(2,1) = -200.0D0*Y(1)*Y(2) - 1.0D0               !
  393. !                                                         !
  394. !      PD(2,2) = 100.0D0*(1.0D0 - Y(1)*Y(1))              !
  395. !                                                         !
  396. !      PD(1,1) = CJ - PD(1,1)                             !
  397. !                                                         !
  398. !      PD(1,2) =    - PD(1,2)                             !
  399. !                                                         !
  400. !      PD(2,1) =    - PD(2,1)                             !
  401. !                                                         !
  402. !      PD(2,2) = CJ - PD(2,2)                             !
  403. !                                                         !
  404. !      RETURN                                             !
  405. !                                                         !
  406. !      END                                                !
  407.  
  408.  
  409.  
  410. gr22=[...
  411. '      SUBROUTINE GR22 (NEQ, T, Y, NG, GROOT, RPAR, IPAR)';
  412. '      IMPLICIT DOUBLE PRECISION (A-H,O-Z)';
  413. '      INTEGER NEQ, NG';
  414. '      DOUBLE PRECISION T, Y, GROOT';
  415. '      DIMENSION Y(*), GROOT(*)';
  416. '      GROOT(1) = Y(1)';
  417. '      RETURN';
  418. '      END';]
  419.  gr22  =
  420.  
  421. !      SUBROUTINE GR22 (NEQ, T, Y, NG, GROOT, RPAR, IPAR)  !
  422. !                                                          !
  423. !      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                 !
  424. !                                                          !
  425. !      INTEGER NEQ, NG                                     !
  426. !                                                          !
  427. !      DOUBLE PRECISION T, Y, GROOT                        !
  428. !                                                          !
  429. !      DIMENSION Y(*), GROOT(*)                            !
  430. !                                                          !
  431. !      GROOT(1) = Y(1)                                     !
  432. !                                                          !
  433. !      RETURN                                              !
  434. !                                                          !
  435. !      END                                                 !
  436.  
  437.  
  438. //Uncomment lines below: link may be machine dependent if some f77 libs are
  439.  
  440. //needed for linking
  441.  
  442. //unix_g('cd /tmp;rm -f /tmp/res22.f');unix_g('cd /tmp;rm -f /tmp/gr22.f');
  443.  
  444. //unix_g('cd /tmp;rm -f /tmp/jac22.f');
  445.  
  446. //write('/tmp/res22.f',res22);write('/tmp/gr22.f',gr22);write('/tmp/jac22.f',jac22)
  447.  
  448. //unix_g("cd /tmp;make /tmp/res22.o");unix_g('cd /tmp;make /tmp/gr22.o');
  449.  
  450. //unix_g('cd /tmp;make /tmp/jac22.o');
  451.  
  452. //          2  Linking the routines
  453.  
  454. //link('/tmp/res22.o','res22');link('/tmp/jac22.o','jac22');link('/tmp/gr22.o','gr22')
  455.  
  456. //rtol=[1.d-6;1.d-6];atol=[1.d-6;1.d-4];
  457.  
  458. //t0=0;y0=[2;0];y0d=[0;-2];t=[20:20:200];ng=1;
  459.  
  460. //info=list([],0,[],[],[],0,0);
  461.  
  462. //          3 Calling the routines by dasrt
  463.  
  464. //[yy,nn]=dasrt([y0,y0d],t0,t,atol,rtol,'res22','jac22',ng,'gr22',info);
  465.  
  466. // hot restart
  467.  
  468. //[yy,nn,hot]=dasrt([y0,y0d],t0,t,atol,rtol,'res22','jac22',ng,'gr22',info);
  469.  
  470. //t01=nn(1);t=100:20:200;[pp,qq]=size(yy);y01=yy(2:3,qq);y0d1=yy(3:4,qq);
  471.  
  472. //[yy,nn,hot]=dasrt([y01,y0d1],t01,t,atol,rtol,'res22','jac22',ng,'gr22',info,hot);
  473.  
  474.  
  475.  
  476.